function filter_ar = get_ar(autocor,p,q)
%from Lesage 2008
%modified by Josh Crozier 2020

% Estimates the AR filter
option = 1;     % Returns all the roots
%option = 2;     % Eliminates the conjugated complex roots

% Takes the middle of the signal
midle=ceil(length(autocor)/2);

% Construction of the autocorrelation matrix

     matautocor=zeros(p,p);

     for k=1:p  
        matautocor(k,:)=autocor(midle-q-k+1:midle-q+p-k);
     end

% Makes the second member vector
secondmember=autocor(midle+q+1:midle+q+p);

% Calculates the coefficients of the associated AR filter
coeffilterAR=matautocor\secondmember';

% Calculates the AR filter poles
coefpolyno=[1 -coeffilterAR'];
polesAR=roots(coefpolyno);
nbroot=size(polesAR,1);

% Check if some poles are out of the unit circle
if q==0
    for k=1:nbroot
         if abs(polesAR(k))>1.6
          error('Some poles are out of the unit circle') 
         end
    end
end %if q

if option == 1
    filter_ar = polesAR;
    
else
    
    realroot=zeros(nbroot,1);
    cplxroot=zeros(nbroot,1);
    nrealroot=0;
    ncplxroot=0;


    % Separates real and complex roots

    for k=1:nbroot

         % Real root
         if (isreal(polesAR(k))==1)
             nrealroot=(nrealroot+1);
             realroot(nrealroot)=polesAR(k);

             % Complex one
         else
             ncplxroot=(ncplxroot+1);
             cplxroot(ncplxroot)=polesAR(k);

         end %endif

    end %endfor

    cplxroot=cplxroot(1:ncplxroot);
    realroot=realroot(1:nrealroot);

    rootcplexsconj=zeros((ncplxroot)/2,1);

    % Eliminates the conjugated complex roots
    for k=1:2:ncplxroot
        rootcplexsconj((k+1)/2)=cplxroot(k);
    end

    filter_ar =[rootcplexsconj' realroot'];

end
end